
library(tidyr)
library(statnet)

###################
### DATA CLEANING
###################
load("Addhealth data with W1 W2 network.RData")


# clean race variable
data$race<-"Black"
data$race[data$latino==1]<-"Other" 
data$race[data$white==1]<-"White"
data$race[data$n.american==1]<-"N. American"
data$race[data$asian==1]<-"Asian"
data$race[data$other==1]<-"Other"

# change so that sex is coded as 0 and 1
data$sex<-data$sex-1

# node level data: 
addhealth_nodes<-data[,c("aid","scid","latino","white","black",
                         "n.american","asian","other","race",
                         "grade","sex", "gpa1", "gpa2")]


##############################
####    CREATE NODESET    ####
##############################
# remove nodes that did not complete the survey in both waves
# include all friendship variables and gpa and grade
# keep GPA as a check because they may not have listed friends, but may have included GPA
node_check <- data[,c("aid", "gpa1", "gpa2","grade", "MF_AID1","MF_AID2","MF_AID3",
                      "MF_AID4","MF_AID5","FF_AID1","FF_AID2","FF_AID3",
                      "FF_AID4","FF_AID5","mf1aid","mf2aid","mf3aid",
                      "mf4aid","mf5aid","ff1aid","ff2aid","ff3aid",
                      "ff4aid","ff5aid")]

var_w1 <- c("gpa1","mf1aid","mf2aid","mf3aid",
            "mf4aid","mf5aid","ff1aid","ff2aid","ff3aid",
            "ff4aid","ff5aid")

var_w2 <- c("gpa2", "MF_AID1","MF_AID2","MF_AID3",
            "MF_AID4","MF_AID5","FF_AID1","FF_AID2","FF_AID3",
            "FF_AID4","FF_AID5")

node_check$missing1 <- ifelse(
  rowSums(is.na(node_check[, var_w1])) == length(var_w1), 1, 0)
# 20 students missing data for all of the variables in wave 1

node_check$missing2 <- ifelse(
  rowSums(is.na(node_check[, var_w2])) == length(var_w2), 1, 0)
# 9 students missing data for all of the variables in wave 2

# nodes to remove: 
nodes_missing <- node_check[(node_check$missing1 == 1 | node_check$missing2 == 1), ]
# these are the 29 nodes to remove - no evidence that they did the survey both waves

nodes_to_remove <- nodes_missing$aid
data_filtered <- data[!(data[[1]] %in% nodes_to_remove), ]
# 1138 rows left

# check for missingness in control variables: 
table(data_filtered$sex, useNA = "always")
# none missing on this variable

table(data_filtered$race, useNA = "always")
# none missing on this variable

table(data_filtered$grade, useNA = "always")
# 22 missing on this variable


# drop students with grade missing
data_filtered <- data_filtered[!is.na(data_filtered$grade),]
nrow(data_filtered)
# 1116 rows left

# drop students with grade <10
data_filtered <- data_filtered[data_filtered$grade>=10,]
# 11 students dropped

# final set of nodes with all controls and having completed the survey in both waves
node_list <- data_filtered[,c("aid", "race", "grade", "sex")]
# 1,105

###############################
####    CREATE EDGELIST    ####
###############################
#### EDGELIST FOR EACH WAVE
# the nodes in the dataset are: 
nodes <- as.integer(data_filtered$aid)

# create edgelist for wave 1:
wave_1<-data_filtered[,c("aid","scid","mf1aid","mf2aid","mf3aid",
                         "mf4aid","mf5aid","ff1aid","ff2aid","ff3aid",
                         "ff4aid","ff5aid")]

el1 <- wave_1 %>%
  pivot_longer(!c('aid', 'scid'), values_to = "alter")


# create flag for AddHealth missing variables - denotes students at other schools
el1$missing_flag <- ifelse(el1$alter %in% c(77777777, 88888888, 99959995, 99999999), el1$alter, ifelse(is.na(el1$alter), "empty", "nonmissing"))

# flag to keep anything that is any version of missing OR is sent to someone in the nodeset
el1$keep <- ifelse(el1$missing_flag == "nonmissing", ifelse(el1$alter %in% nodes, 1, 0), 1)
# get rid of edges that are sent to students who are not in the dataset
el1_clean <- el1[el1$keep == 1, ]

# filter for only actual edges:
el1_nonmissing <- el1[(el1$keep == 1 & el1$missing_flag == "nonmissing"), ]

missing_nodes_w1 <- c(77777777, 88888888, 99959995, 99999999)

# now, go back to el1_clean, which has the edgelist with edges sent to students who we don't have data for removed
# what is remaining here: 
## edges between students for whom we have complete data
## edges from a student with complete data to one coded as out of school
## edges from a student for whom we have complete data to NA
total_nodes_w1 = c(as.numeric(el1_clean$aid), el1_clean$alter)
total_w1 = unique(total_nodes_w1)
total_w1 = total_w1[!(total_w1 %in% missing_nodes_w1)]
total_w1 = total_w1[!is.na(total_w1)]
length(total_w1)
# 1,105 node in the network with the lower grades removed


# create edgelist for wave 2:
wave_2<-data_filtered[,c("aid","scid","MF_AID1","MF_AID2","MF_AID3",
                         "MF_AID4","MF_AID5","FF_AID1","FF_AID2","FF_AID3",
                         "FF_AID4","FF_AID5")]

# convert to edgelist for Wave 2
el2 <- wave_2 %>%
  pivot_longer(!c('aid', 'scid'), values_to = "alter")

# create flag for AddHealth missing variables - denotes students at other schools
el2$missing_flag <- ifelse(el2$alter %in% c(77777777, 88888888, 99959995, 99999999, 55555555), el2$alter, ifelse(is.na(el2$alter), "empty", 'nonmissing'))
# flag to keep anything that is missing OR is sent to someone in the nodeset
el2$keep <- ifelse(el2$missing_flag == "nonmissing", ifelse(el2$alter %in% nodes, 1, 0), 1)
# get rid of edges that are sent to students who have been removed from dataset
el2_clean <- el2[el2$keep == 1, ]
# only kept edges that are missing or within the nodeset

# if we filter for only actual edges, we get:
el2_nonmissing <- el2[(el2$keep == 1 & el2$missing_flag == "nonmissing"), ]

# now, see how many nodes exist in network:
missing_nodes_w2 <- c(77777777, 88888888, 99959995, 99999999, 55555555)
total_nodes_w2 = c(as.numeric(el2_clean$aid), el2_clean$alter)
total_w2 = unique(total_nodes_w2)
total_w2 = total_w2[!(total_w2 %in% missing_nodes_w2)]
total_w2 = total_w2[!is.na(total_w2)]
length(total_w2)
# 1,105 nodes in this network


###############################################
## REMOVE LOOPS AND DUPLICATES FROM EDGELISTS
###############################################
# wave 1
el1_nonmissing
wave_1_el <- el1_nonmissing[,c("aid", "alter")]

duplicates <- wave_1_el[duplicated(wave_1_el), ]
# no duplicates here

# also, need to remove loops
# check for and remove loops:
wave_1_el$check <- ifelse(wave_1_el$aid == wave_1_el$alter, 1, 0)
wave_1_el <- wave_1_el[wave_1_el$check == 0,]

# now, get rid of the check column:
wave_1_el <- wave_1_el[,c("aid", "alter")]



# wave 2
el2_nonmissing
wave_2_el <- el2_nonmissing[,c("aid", "alter")]
duplicates <- wave_2_el[duplicated(wave_2_el), ]

# drop duplicates
wave_2_el <- wave_2_el[!duplicated(wave_2_el), ]

# also, need to remove loops
# check for and remove loops:
wave_2_el$check <- ifelse(wave_2_el$aid == wave_2_el$alter, 1, 0)
wave_2_el <- wave_2_el[wave_2_el$check == 0,]

# now, get rid of the check column:
wave_2_el <- wave_2_el[,c("aid", "alter")]



##################
## BUILD NETWORK
##################
# inputs to network
# wave 1 edgelist: wave_1_el
# wave 2 edgelist: wave_2_el
# nodes: node_list

# first, create a network object from these inputs:
library(igraph)
library(intergraph)
library(statnet)

g_1 = graph_from_data_frame(d = wave_1_el, vertices = node_list)
g_2 = graph_from_data_frame(d = wave_2_el, vertices = node_list)


net_1 <- asNetwork(g_1)
# density is 0.001645049

# check percent of non-null dyads that are reciprocated:
grecip(net_1, measure = "dyadic.nonnull")
# 0.2058997 

net_2 <- asNetwork(g_2)
# density is 0.001305913

# check percent of non-null dyads that are reciprocated:
grecip(net_2, measure = "dyadic.nonnull")
# 0.2065623 

# check global clustering: 
gtrans(net_1, mode = "global")
# 0.2387199
gtrans(net_2, mode = "global")
# 0.236799


#####################
### TERGM
#####################
# put networks into a list
net_list <- list()
net_list[[1]] <- net_1
net_list[[2]] <- net_2


install.packages("ergm.multi", type = "source")
library(tergm)
TERGM_1<-tergm(net_list~Form(~edges + 
                               mutual + 
                               nodeifactor("sex") + 
                               nodeofactor("sex") + 
                               nodematch("sex") + 
                               nodeifactor("race", levels = -1) + 
                               nodeofactor("race", levels = -1) + 
                               nodematch("race", diff=TRUE, levels = -1) + 
                               absdiff("grade")),
                               estimate="CMLE",
                               control=control.tergm(parallel=3) #run on 3 processors to speed up estimation
)

summary(TERGM_1)


# add in transitivity term:
TERGM_2<-tergm(net_list~Form(~edges + 
                               mutual + 
                               nodeifactor("sex") + 
                               nodeofactor("sex") + 
                               nodematch("sex") + 
                               nodeifactor("race", levels = -1) + 
                               nodeofactor("race", levels = -1) + 
                               nodematch("race", diff=TRUE, levels = -1) + 
                               absdiff("grade") + 
                               gwesp(.7,fixed=T)),
                               estimate="CMLE",
                               control=control.tergm(parallel=3) #run on 3 processors to speed up estimation
)

summary(TERGM_2)


######################
###  GOF - TERGM  ####
######################
library(ergMargins)

# read in models
TERGM_1 <- readRDS("TERGM_1_21_04__11_23.rds")
TERGM_2 <- readRDS("TERGM_2_21_04__12_11.rds")

# note: TERGM_1 didn't have endogenous term so it did not use MCMC
## check MCMC mixing
mcmc.diagnostics(TERGM_2)
par(mar=c(1,1,1,1))
plot(TERGM_2$sample) ##plot MCMC diagnostics


##check goodness of fit
gof.TERGM_1<-gof(TERGM_1) 
gof.TERGM_1$pval.model
gof.TERGM_1$summary.model

ergm:::plot.gof(gof.TERGM_1,roc.col="black",pr.col="black",obs.col="black") 


gof.TERGM_2<-gof(TERGM_2) 
plot(gof.TERGM_2,roc.col="black",pr.col="black",obs.col="black")

#check collinearity
ergMargins::vif.ergm(TERGM_1)  ##calculate variance inflation factors
ergMargins::vif.ergm(TERGM_2)  ##calculate variance inflation factors


#################################
##### SIENA MODELS      #########
#################################
install.packages("RSiena") 
library(RSiena)

# first, get adjacency matrix from the statnet objects
adj_net_1 <- as.matrix(net_1, matrix.type = "adjacency")
adj_net_2 <- as.matrix(net_2, matrix.type = "adjacency")


##create "sienaDependent" object
AddHealthNet<-sienaDependent(array(c(adj_net_1,
                                     adj_net_2),
                                   dim=c(1105,1105,2)))

###Now assign variables
##get sex attribute
sex<-net_1%v%"sex"
##sex is non-time varying, so we set the attribute using "coCovar"
Sex<-coCovar(sex)

#first, convert race to three binary variables - omitted category is Asian
white <- ifelse(net_1%v%"race" == "White", 1, 0)
black <- ifelse(net_1%v%"race" == "Black", 1, 0)
other <- ifelse(net_1%v%"race" == "Other", 1, 0)

# now, get each attribute and set it as a coCovar
White<-coCovar(white)
Black<-coCovar(black)
Other<-coCovar(other)

#get grade attribute
grade<-net_1%v%"grade"
Grade<-coCovar(grade)

#create the SIENA dataset
SAOM.Data<-sienaDataCreate(Network=AddHealthNet,Sex,White,Black,Other,Grade)


#  formulate the model
SAOM.terms<-getEffects(SAOM.Data)    

SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,sameX,interaction1="White") 
SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,sameX,interaction1="Black") 
SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,sameX,interaction1="Other") 
SAOM.terms<-includeEffects(SAOM.terms,egoX,altX,sameX,interaction1="Sex") 
SAOM.terms<-includeEffects(SAOM.terms,absDiffX,interaction1="Grade") 

# create the model
create.model<-sienaAlgorithmCreate(projname="AddHealth_1", 
                                   seed=21093, 
                                   nsub=5,   
                                   n3=2000) 


# estimate the model using siena07
AddHealthModel_1<-siena07(create.model,
                          data=SAOM.Data,
                          effects=SAOM.terms,
                          verbose=T,  
                          returnDeps=T) 

AddHealthModel_1


# for the next model, include transitive triplets:
SAOM.terms_2<-includeEffects(SAOM.terms,transTies)

# create the model
create.model_2<-sienaAlgorithmCreate(projname="AddHealth_2", 
                                     seed=21093, 
                                     nsub=5,  
                                     n3=2000) 


# estimate the model using siena07
AddHealthModel_2<-siena07(create.model_2,
                          data=SAOM.Data,
                          effects=SAOM.terms_2,
                          verbose=T,  
                          returnDeps=T) 

AddHealthModel_2


######################
###  GOF - SIENA  ####
######################
# AddHealth p-values 
# Model 1
Siena_1_results <- data.frame(var = AddHealthModel_1$effects$effectName, theta = AddHealthModel_1$theta, se = AddHealthModel_1$se)
Siena_1_results$z <- Siena_1_results$theta / Siena_1_results$se
Siena_1_results$p_value <- 2 * (1 - pnorm(abs(Siena_1_results$z)))

# Model 2
Siena_2_results <- data.frame(var = AddHealthModel_2$effects$effectName, theta = AddHealthModel_2$theta, se = AddHealthModel_2$se)
Siena_2_results$z <- Siena_2_results$theta / Siena_2_results$se
Siena_2_results$p_value <- 2 * (1 - pnorm(abs(Siena_2_results$z)))



###Assess goodness of fit
indegGOF<-sienaGOF(AddHealthModel_1,IndegreeDistribution,join = T,varName="Network")
outdegGOF<-sienaGOF(AddHealthModel_1,OutdegreeDistribution,join = T,varName="Network")
TriadGOF<-sienaGOF(AddHealthModel_1,TriadCensus,join = T,varName="Network")


plot(indegGOF)
plot(outdegGOF)
plot(TriadGOF)


indegGOF<-sienaGOF(AddHealthModel_2,IndegreeDistribution,join = T,varName="Network")
outdegGOF<-sienaGOF(AddHealthModel_2,OutdegreeDistribution,join = T,varName="Network")
TriadGOF<-sienaGOF(AddHealthModel_2,TriadCensus,join = T,varName="Network")


plot(indegGOF)
plot(outdegGOF)
plot(TriadGOF)


############################
#####  MEMS - ALL   #####
############################
### MACRO STRUCTURE: COLEMAN SEGREGATION INDEX
# the function is written directly into the MEMS calculations, but need group size for all calculations
# create group_size variable: 
race_vector <- network::get.vertex.attribute(net_1, "race")
asian_group = table(race_vector)[[1]]
black_group = table(race_vector)[[2]]
other_group = table(race_vector)[[3]]
white_group = table(race_vector)[[4]]
group_size <- c(asian_group, black_group, other_group, white_group)


# MEMS for TERGM_1 - in-group preference
MEMS_1 = MEMS(TERGM_1,
              micro_process = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other"),
              macro_function = function(net){
                mixing_matrix <- mixingmatrix(net, "race")
                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                return(coleman_index[1])
              } ,
              object_type = 'network',
              silent=F,
              algorithm ="parametric")


# MEMS for TERGM_1 - avoidance by other students:
MEMS_2 = MEMS(TERGM_1,
              micro_process = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other", "Form(1)~nodematch.race.Black", "Form(1)~nodematch.race.Other", "Form(1)~nodematch.race.White"),
              macro_function = function(net){
                mixing_matrix <- mixingmatrix(net, "race")
                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                return(coleman_index[1])
              } ,
              object_type = 'network',
              silent=F,
              algorithm ="parametric")


# compare TERGM MEMS for in-group preference:
change_MEMS_1 <- compare_MEMS(partial_model = TERGM_1, 
                              full_model = TERGM_2,
                              micro_process = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other"),
                              macro_function = function(net){
                                mixing_matrix <- mixingmatrix(net, "race")
                                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                return(coleman_index[1])
                              },
                              silent=F,
                              algorithm ="parametric")



# compare TERGM MEMS for avoidance by other students:
change_MEMS_2 <- compare_MEMS(partial_model = TERGM_1, 
                              full_model = TERGM_2, 
                              micro_process = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other", "Form(1)~nodematch.race.Black", "Form(1)~nodematch.race.Other", "Form(1)~nodematch.race.White"),
                              macro_function =function(net){
                                mixing_matrix <- mixingmatrix(net, "race")
                                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                return(coleman_index[1])
                              },
                              silent=F,
                              algorithm ="parametric")


# MEMS for Siena_1 - in-group preference
MEMS_3 = MEMS(AddHealthModel_1,
              micro_process = c("Black alter","White alter", "Other alter"),
              macro_function = function(net){
                if(!"race"%in%network::list.vertex.attributes(net)){
                  race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                  race_df[race_df > 0] <- 1 
                  race_df[race_df < 0] <- 0
                  race_df$race <- "Asian"
                  race_df$race[race_df$black == 1] <- "Black"
                  race_df$race[race_df$white == 1] <- "White"
                  race_df$race[race_df$other == 1] <- "Other"
                  
                  network::set.vertex.attribute(net,"race",race_df$race)
                  mixing_matrix <- mixingmatrix(net, "race")
                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                  return(coleman_index[1])
                } else{
                  mixing_matrix <- mixingmatrix(net, "race")
                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                  return(coleman_index[1])
                }
              },
              nsim=500,
              SAOM_data = SAOM.Data,
              object_type = 'network',
              silent=F,
              algorithm ="parametric")



# MEMS for Siena_1 - avoidance by other students:
MEMS_4 = MEMS(AddHealthModel_1,
              micro_process = c("Black alter","White alter", "Other alter", "same Black", "same White", "same Other"),
              macro_function = function(net){
                # first, I need to create a race attribute based on the individual race labels from the network:
                # create a race vector:
                # race vector:
                if(!"race"%in%network::list.vertex.attributes(net)){
                  race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                  race_df[race_df > 0] <- 1 
                  race_df[race_df < 0] <- 0
                  race_df$race <- "Asian"
                  race_df$race[race_df$black == 1] <- "Black"
                  race_df$race[race_df$white == 1] <- "White"
                  race_df$race[race_df$other == 1] <- "Other"
                  
                  network::set.vertex.attribute(net,"race",race_df$race)
                  mixing_matrix <- mixingmatrix(net, "race")
                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                  return(coleman_index[1])
                } else{
                  mixing_matrix <- mixingmatrix(net, "race")
                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                  return(coleman_index[1])
                }
              },
              nsim=500,
              SAOM_data = SAOM.Data,
              object_type = 'network',
              silent=F,
              algorithm ="parametric")



# compare Siena MEMS for in-group preference:
change_MEMS_3 <- compare_MEMS(partial_model = AddHealthModel_1, 
                              full_model = AddHealthModel_2, 
                              micro_process = c("Black alter","White alter", "Other alter"),
                              macro_function =function(net){
                                # first, I need to create a race attribute based on the individual race labels from the network:
                                # create a race vector:
                                # race vector:
                                if(!"race"%in%network::list.vertex.attributes(net)){
                                  race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                  race_df[race_df > 0] <- 1 
                                  race_df[race_df < 0] <- 0
                                  race_df$race <- "Asian"
                                  race_df$race[race_df$black == 1] <- "Black"
                                  race_df$race[race_df$white == 1] <- "White"
                                  race_df$race[race_df$other == 1] <- "Other"
                                  
                                  network::set.vertex.attribute(net,"race",race_df$race)
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                } else{
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                }
                              },
                              full_output=T,
                              SAOM_data = SAOM.Data,
                              silent=F,
                              algorithm ="parametric")

# compare Siena MEMS for out-group avoidance:
change_MEMS_4 <- compare_MEMS(partial_model = AddHealthModel_1, 
                              full_model = AddHealthModel_2, 
                              micro_process = c("Black alter","White alter", "Other alter", "same Black", "same White", "same Other"),
                              macro_function =function(net){
                                if(!"race"%in%network::list.vertex.attributes(net)){
                                  race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                  race_df[race_df > 0] <- 1 
                                  race_df[race_df < 0] <- 0
                                  race_df$race <- "Asian"
                                  race_df$race[race_df$black == 1] <- "Black"
                                  race_df$race[race_df$white == 1] <- "White"
                                  race_df$race[race_df$other == 1] <- "Other"
                                  
                                  network::set.vertex.attribute(net,"race",race_df$race)
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                } else{
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                }
                              },
                              full_output=T,
                              SAOM_data = SAOM.Data,
                              silent=F,
                              algorithm ="parametric")


# sensitivity test for functional form:
# first, for in-group preference:
change_MEMS_5 <- compare_MEMS(full_model = TERGM_1, 
                              partial_model = AddHealthModel_1, 
                              micro_process2 = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other"),
                              micro_process = c("Black alter","White alter", "Other alter"),
                              macro_function =function(net){
                                if(!"race"%in%network::list.vertex.attributes(net)){
                                  race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                  race_df[race_df > 0] <- 1 
                                  race_df[race_df < 0] <- 0
                                  race_df$race <- "Asian"
                                  race_df$race[race_df$black == 1] <- "Black"
                                  race_df$race[race_df$white == 1] <- "White"
                                  race_df$race[race_df$other == 1] <- "Other"
                                  
                                  network::set.vertex.attribute(net,"race",race_df$race)
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                } else{
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                }
                              },
                              full_output=T,
                              SAOM_data = SAOM.Data,
                              silent=F,
                              algorithm ="parametric")

# next, for in-group preference and inclusion of transitive triads
change_MEMS_6 <- compare_MEMS(full_model = TERGM_2, 
                              partial_model = AddHealthModel_2, 
                              micro_process2 = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other"),
                              micro_process = c("Black alter","White alter", "Other alter"),
                              macro_function =function(net){
                                if(!"race"%in%network::list.vertex.attributes(net)){
                                  race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                  race_df[race_df > 0] <- 1 
                                  race_df[race_df < 0] <- 0
                                  race_df$race <- "Asian"
                                  race_df$race[race_df$black == 1] <- "Black"
                                  race_df$race[race_df$white == 1] <- "White"
                                  race_df$race[race_df$other == 1] <- "Other"
                                  
                                  network::set.vertex.attribute(net,"race",race_df$race)
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                } else{
                                  mixing_matrix <- mixingmatrix(net, "race")
                                  coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                  return(coleman_index[1])
                                }
                              },
                              full_output=T,
                              SAOM_data = SAOM.Data,
                              silent=F,
                              algorithm ="parametric")

# now for out-group avoidance, no transitive triads
change_MEMS_7<- compare_MEMS(full_model = TERGM_1, 
                             partial_model = AddHealthModel_1, 
                             micro_process2 = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other", "Form(1)~nodematch.race.Black", "Form(1)~nodematch.race.Other", "Form(1)~nodematch.race.White"),
                             micro_process = c("Black alter","White alter", "Other alter", "same Black", "same White", "same Other"),
                             macro_function =function(net){
                               if(!"race"%in%network::list.vertex.attributes(net)){
                                 race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                 race_df[race_df > 0] <- 1 
                                 race_df[race_df < 0] <- 0
                                 race_df$race <- "Asian"
                                 race_df$race[race_df$black == 1] <- "Black"
                                 race_df$race[race_df$white == 1] <- "White"
                                 race_df$race[race_df$other == 1] <- "Other"
                                 
                                 network::set.vertex.attribute(net,"race",race_df$race)
                                 mixing_matrix <- mixingmatrix(net, "race")
                                 coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                 return(coleman_index[1])
                               } else{
                                 mixing_matrix <- mixingmatrix(net, "race")
                                 coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                 return(coleman_index[1])
                               }
                             },
                             full_output=T,
                             SAOM_data = SAOM.Data,
                             silent=F,
                             algorithm ="parametric")

# out-group avoidance, inclusion of transitive triads
change_MEMS_8<- compare_MEMS(full_model = TERGM_2, 
                             partial_model = AddHealthModel_2, 
                             micro_process2 = c("Form(1)~nodeifactor.race.Black","Form(1)~nodeifactor.race.White", "Form(1)~nodeifactor.race.Other", "Form(1)~nodematch.race.Black", "Form(1)~nodematch.race.Other", "Form(1)~nodematch.race.White"),
                             micro_process = c("Black alter","White alter", "Other alter", "same Black", "same White", "same Other"),
                             macro_function =function(net){
                               if(!"race"%in%network::list.vertex.attributes(net)){
                                 race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                 race_df[race_df > 0] <- 1 
                                 race_df[race_df < 0] <- 0
                                 race_df$race <- "Asian"
                                 race_df$race[race_df$black == 1] <- "Black"
                                 race_df$race[race_df$white == 1] <- "White"
                                 race_df$race[race_df$other == 1] <- "Other"
                                 
                                 network::set.vertex.attribute(net,"race",race_df$race)
                                 mixing_matrix <- mixingmatrix(net, "race")
                                 coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                 return(coleman_index[1])
                               } else{
                                 mixing_matrix <- mixingmatrix(net, "race")
                                 coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                 return(coleman_index[1])
                               }
                             },
                             full_output=T,
                             SAOM_data = SAOM.Data,
                             silent=F,
                             algorithm ="parametric")


### MEMS for Siena model:
# MEMS for Siena without triadic closure
AddHealthModel_1_MEMS <- list()
for(i in 1:length(AddHealthModel_1$effects$effectName)){
  AddHealthModel_1_MEMS[[i]] <- MEMS(AddHealthModel_1,
                                     micro_process = AddHealthModel_1$effects$effectName[i],
                                     macro_function = function(net){
                                       # first, I need to create a race attribute based on the individual race labels from the network:
                                       # create a race vector:
                                       # race vector:
                                       if(!"race"%in%network::list.vertex.attributes(net)){
                                         race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                         race_df[race_df > 0] <- 1 
                                         race_df[race_df < 0] <- 0
                                         race_df$race <- "Asian"
                                         race_df$race[race_df$black == 1] <- "Black"
                                         race_df$race[race_df$white == 1] <- "White"
                                         race_df$race[race_df$other == 1] <- "Other"
                                         
                                         network::set.vertex.attribute(net,"race",race_df$race)
                                         mixing_matrix <- mixingmatrix(net, "race")
                                         coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                         return(coleman_index[1])
                                       } else{
                                         mixing_matrix <- mixingmatrix(net, "race")
                                         coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                         return(coleman_index[1])
                                       }
                                     },
                                     nsim=500,
                                     SAOM_data = SAOM.Data,
                                     object_type = 'network',
                                     silent=F,
                                     algorithm ="parametric")
}



# MEMS for Siena 1 - avoidance by other students:
AddHealthModel_2_MEMS <- list()
for(i in 1:length(AddHealthModel_2$effects$effectName)){
  AddHealthModel_2_MEMS[[i]] <- MEMS(AddHealthModel_2,
                                     micro_process = AddHealthModel_2$effects$effectName[i],
                                     macro_function = function(net){
                                       # first, I need to create a race attribute based on the individual race labels from the network:
                                       # create a race vector:
                                       # race vector:
                                       if(!"race"%in%network::list.vertex.attributes(net)){
                                         race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                         race_df[race_df > 0] <- 1 
                                         race_df[race_df < 0] <- 0
                                         race_df$race <- "Asian"
                                         race_df$race[race_df$black == 1] <- "Black"
                                         race_df$race[race_df$white == 1] <- "White"
                                         race_df$race[race_df$other == 1] <- "Other"
                                         
                                         network::set.vertex.attribute(net,"race",race_df$race)
                                         mixing_matrix <- mixingmatrix(net, "race")
                                         coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                         return(coleman_index[1])
                                       } else{
                                         mixing_matrix <- mixingmatrix(net, "race")
                                         coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                         return(coleman_index[1])
                                       }
                                     },
                                     nsim=500,
                                     SAOM_data = SAOM.Data,
                                     object_type = 'network',
                                     silent=F,
                                     algorithm ="parametric")
}




### MEMS TERGM:
# MEMS for TERGM without triadic closure
TERGM_1_MEMS <- list()
for(i in 1:length(names(TERGM_1$coefficients))){
  TERGM_1_MEMS[[i]] <- MEMS(TERGM_1,
                            micro_process = names(TERGM_1$coefficients)[i],
                            macro_function = function(net){
                              if(!"race"%in%network::list.vertex.attributes(net)){
                                race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                race_df[race_df > 0] <- 1 
                                race_df[race_df < 0] <- 0
                                race_df$race <- "Asian"
                                race_df$race[race_df$black == 1] <- "Black"
                                race_df$race[race_df$white == 1] <- "White"
                                race_df$race[race_df$other == 1] <- "Other"
                                
                                network::set.vertex.attribute(net,"race",race_df$race)
                                mixing_matrix <- mixingmatrix(net, "race")
                                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                return(coleman_index[1])
                              } else{
                                mixing_matrix <- mixingmatrix(net, "race")
                                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                return(coleman_index[1])
                              }
                            },
                            nsim=500,
                            SAOM_data = SAOM.Data,
                            object_type = 'network',
                            silent=F,
                            algorithm ="parametric")
}



# MEMS for TERGM without triadic closure
TERGM_2_MEMS <- list()
for(i in 1:length(names(TERGM_2$coefficients))){
  TERGM_2_MEMS[[i]] <- MEMS(TERGM_2,
                            micro_process = names(TERGM_2$coefficients)[i],
                            macro_function = function(net){
                              if(!"race"%in%network::list.vertex.attributes(net)){
                                race_df <- data.frame(black = c(network::get.vertex.attribute(net, "Black")), white = c(network::get.vertex.attribute(net, "White")), other = c(network::get.vertex.attribute(net, "Other")))
                                race_df[race_df > 0] <- 1 
                                race_df[race_df < 0] <- 0
                                race_df$race <- "Asian"
                                race_df$race[race_df$black == 1] <- "Black"
                                race_df$race[race_df$white == 1] <- "White"
                                race_df$race[race_df$other == 1] <- "Other"
                                
                                network::set.vertex.attribute(net,"race",race_df$race)
                                mixing_matrix <- mixingmatrix(net, "race")
                                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                return(coleman_index[1])
                              } else{
                                mixing_matrix <- mixingmatrix(net, "race")
                                coleman_index <- coleman(mixing_matrix, gsizes = group_size)
                                return(coleman_index[1])
                              }
                            },
                            nsim=500,
                            SAOM_data = SAOM.Data,
                            object_type = 'network',
                            silent=F,
                            algorithm ="parametric")
}


###################################################
#### compare the MEMS between models
###################################################
# Take the samples created in the MEMS calculation and comparing them across models

# rename models:
in_group_preference <- MEMS_3
out_group_avoidance <- MEMS_4
change_MEMS_in_group <- change_MEMS_3
change_MEMS_out_group <- change_MEMS_4
change_MEMS_in_group_TERGM <- change_MEMS_1
change_MEMS_out_group_TERGM <- change_MEMS_2
change_MEMS_in_group_specification <- change_MEMS_5
change_MEMS_out_group_specification <- change_MEMS_7
change_MEMS_in_group_specification_triadic <- rchange_MEMS_6
change_MEMS_out_group_specification_triadic <- change_MEMS_8


###########################
### IN-GROUP AVOIDANCE ###
###########################
# create nsim variable from the MEMS calculation
nsim=500

## INDIRECT EFFECT
# this is the change in MEMS for the Siena model
change_MEMS_in_group$diff_MEMS_results$diff_samples

# this is the change in MEMS for the TERGM
change_MEMS_in_group_TERGM$diff_MEMS_results$diff_samples

diff_diff_MEMS_dat<-matrix(NA,nrow=3,ncol=7)
colnames(diff_diff_MEMS_dat)<-c(colnames(change_MEMS_in_group$diff_MEMS_results$diff_summary),"Prop. of indirect MEMS explained")

r_names<-c("Indirect MEMS - SAOM","Indirect MEMS - TERGM","Diff. in Indirect MEMS")
rownames(diff_diff_MEMS_dat)<-r_names

diff_diff_MEMS_dat[1,]<-c(change_MEMS_in_group$diff_MEMS_results$diff_summary[3,],NA)
diff_diff_MEMS_dat[2,]<-c(change_MEMS_in_group_TERGM$diff_MEMS_results$diff_summary[3,],NA)

diff_diff_MEMS <- change_MEMS_in_group$diff_MEMS_results$diff_samples - change_MEMS_in_group_TERGM$diff_MEMS_results$diff_samples
diff_diff_MEMS_dat[3,1]<-mean(diff_diff_MEMS)
diff_diff_MEMS_dat[3,2]<-sd(diff_diff_MEMS)
diff_diff_MEMS_dat[3,3]<-quantile(diff_diff_MEMS,.025,na.rm=TRUE)
diff_diff_MEMS_dat[3,4]<-quantile(diff_diff_MEMS,.975,na.rm=TRUE)
if(diff_diff_MEMS_dat[3,1]<0){
  diff_diff_MEMS_dat[3,5]<-length(diff_diff_MEMS[which(diff_diff_MEMS>0)])/(nsim)
}else{
  diff_diff_MEMS_dat[3,5]<-length(diff_diff_MEMS[which(diff_diff_MEMS<0)])/(nsim)
  
}
diff_diff_MEMS_dat[3,7]<-mean(diff_diff_MEMS)/change_MEMS_in_group$diff_MEMS_results$diff_summary[3,1]

indirect_MEMS_in <- diff_diff_MEMS_dat


## TOTAL EFFECT
# this is the change in MEMS for the Siena model
change_MEMS_in_group$diff_MEMS_results$diff_samples

# this is the change in MEMS for the TERGM
change_MEMS_in_group_TERGM$diff_MEMS_results$diff_samples

diff_total_MEMS_dat<-matrix(NA,nrow=3,ncol=6)
colnames(diff_total_MEMS_dat)<-c(colnames(change_MEMS_in_group$p_MEMS_results$summary_dat),"Prop. of total MEMS explained")

r_names<-c("Total MEMS - SAOM","Total MEMS - TERGM","Diff. in Total MEMS")
rownames(diff_total_MEMS_dat)<-r_names

diff_total_MEMS_dat[1,]<-c(change_MEMS_in_group$p_MEMS_results$summary_dat[1,],NA)
diff_total_MEMS_dat[2,]<-c(change_MEMS_in_group_TERGM$p_MEMS_results$summary_dat[1,],NA)

diff_total_MEMS <- change_MEMS_in_group$p_MEMS_results$mems_samples - change_MEMS_in_group_TERGM$p_MEMS_results$mems_samples
diff_total_MEMS_dat[3,1]<-mean(diff_total_MEMS)
diff_total_MEMS_dat[3,2]<-sd(diff_total_MEMS)
diff_total_MEMS_dat[3,3]<-quantile(diff_total_MEMS,.025,na.rm=TRUE)
diff_total_MEMS_dat[3,4]<-quantile(diff_total_MEMS,.975,na.rm=TRUE)
if(diff_total_MEMS_dat[3,1]<0){
  diff_total_MEMS_dat[3,5]<-length(diff_total_MEMS[which(diff_total_MEMS>0)])/(nsim)
}else{
  diff_total_MEMS_dat[3,5]<-length(diff_total_MEMS[which(diff_total_MEMS<0)])/(nsim)
  
}
diff_total_MEMS_dat[3,6]<-mean(diff_total_MEMS)/change_MEMS_in_group$p_MEMS_results$summary_dat[1,1]

total_MEMS_in <- diff_total_MEMS_dat


## DIRECT EFFECT

diff_direct_MEMS_dat<-matrix(NA,nrow=3,ncol=6)
colnames(diff_direct_MEMS_dat)<-c(colnames(change_MEMS_in_group$f_MEMS_results$summary_dat),"Prop. of total MEMS explained")

r_names<-c("Direct MEMS - SAOM","Direct MEMS - TERGM","Diff. in Direct MEMS")
rownames(diff_direct_MEMS_dat)<-r_names

diff_direct_MEMS_dat[1,]<-c(change_MEMS_in_group$f_MEMS_results$summary_dat[1,],NA)
diff_direct_MEMS_dat[2,]<-c(change_MEMS_in_group_TERGM$f_MEMS_results$summary_dat[1,],NA)

diff_direct_MEMS <- change_MEMS_in_group$f_MEMS_results$mems_samples - change_MEMS_in_group_TERGM$f_MEMS_results$mems_samples
diff_direct_MEMS_dat[3,1]<-mean(diff_direct_MEMS)
diff_direct_MEMS_dat[3,2]<-sd(diff_direct_MEMS)
diff_direct_MEMS_dat[3,3]<-quantile(diff_direct_MEMS,.025,na.rm=TRUE)
diff_direct_MEMS_dat[3,4]<-quantile(diff_direct_MEMS,.975,na.rm=TRUE)
if(diff_direct_MEMS_dat[3,1]<0){
  diff_direct_MEMS_dat[3,5]<-length(diff_direct_MEMS[which(diff_direct_MEMS>0)])/(nsim)
}else{
  diff_direct_MEMS_dat[3,5]<-length(diff_direct_MEMS[which(diff_direct_MEMS<0)])/(nsim)
  
}
diff_direct_MEMS_dat[3,6]<-mean(diff_direct_MEMS)/change_MEMS_in_group$f_MEMS_results$summary_dat[1,1]

direct_MEMS_in <- diff_direct_MEMS_dat



###########################
### OUT-GROUP AVOIDANCE ###
###########################
## INDIRECT EFFECT
# this is the change in MEMS for the Siena model
change_MEMS_out_group$diff_MEMS_results$diff_samples

# this is the change in MEMS for the TERGM
change_MEMS_out_group_TERGM$diff_MEMS_results$diff_samples

diff_diff_MEMS_dat<-matrix(NA,nrow=3,ncol=7)
colnames(diff_diff_MEMS_dat)<-c(colnames(change_MEMS_out_group$diff_MEMS_results$diff_summary),"Prop. of indirect MEMS explained")

r_names<-c("Indirect MEMS - SAOM","Indirect MEMS - TERGM","Diff. in Indirect MEMS")
rownames(diff_diff_MEMS_dat)<-r_names

diff_diff_MEMS_dat[1,]<-c(change_MEMS_out_group$diff_MEMS_results$diff_summary[3,],NA)
diff_diff_MEMS_dat[2,]<-c(change_MEMS_out_group_TERGM$diff_MEMS_results$diff_summary[3,],NA)

diff_diff_MEMS <- change_MEMS_out_group$diff_MEMS_results$diff_samples - change_MEMS_out_group_TERGM$diff_MEMS_results$diff_samples
diff_diff_MEMS_dat[3,1]<-mean(diff_diff_MEMS)
diff_diff_MEMS_dat[3,2]<-sd(diff_diff_MEMS)
diff_diff_MEMS_dat[3,3]<-quantile(diff_diff_MEMS,.025,na.rm=TRUE)
diff_diff_MEMS_dat[3,4]<-quantile(diff_diff_MEMS,.975,na.rm=TRUE)
if(diff_diff_MEMS_dat[3,1]<0){
  diff_diff_MEMS_dat[3,5]<-length(diff_diff_MEMS[which(diff_diff_MEMS>0)])/(nsim)
}else{
  diff_diff_MEMS_dat[3,5]<-length(diff_diff_MEMS[which(diff_diff_MEMS<0)])/(nsim)
  
}
diff_diff_MEMS_dat[3,7]<-mean(diff_diff_MEMS)/change_MEMS_out_group$diff_MEMS_results$diff_summary[3,1]

indirect_MEMS_out <- diff_diff_MEMS_dat



## TOTAL EFFECT
# this is the change in MEMS for the Siena model
change_MEMS_out_group$diff_MEMS_results$diff_samples

# this is the change in MEMS for the TERGM
change_MEMS_out_group_TERGM$diff_MEMS_results$diff_samples

diff_total_MEMS_dat<-matrix(NA,nrow=3,ncol=6)
colnames(diff_total_MEMS_dat)<-c(colnames(change_MEMS_out_group$p_MEMS_results$summary_dat),"Prop. of total MEMS explained")

r_names<-c("Total MEMS - SAOM","Total MEMS - TERGM","Diff. in Total MEMS")
rownames(diff_total_MEMS_dat)<-r_names

diff_total_MEMS_dat[1,]<-c(change_MEMS_out_group$p_MEMS_results$summary_dat[1,],NA)
diff_total_MEMS_dat[2,]<-c(change_MEMS_out_group_TERGM$p_MEMS_results$summary_dat[1,],NA)

diff_total_MEMS <- change_MEMS_out_group$p_MEMS_results$mems_samples - change_MEMS_out_group_TERGM$p_MEMS_results$mems_samples
diff_total_MEMS_dat[3,1]<-mean(diff_total_MEMS)
diff_total_MEMS_dat[3,2]<-sd(diff_total_MEMS)
diff_total_MEMS_dat[3,3]<-quantile(diff_total_MEMS,.025,na.rm=TRUE)
diff_total_MEMS_dat[3,4]<-quantile(diff_total_MEMS,.975,na.rm=TRUE)
if(diff_total_MEMS_dat[3,1]<0){
  diff_total_MEMS_dat[3,5]<-length(diff_total_MEMS[which(diff_total_MEMS>0)])/(nsim)
}else{
  diff_total_MEMS_dat[3,5]<-length(diff_total_MEMS[which(diff_total_MEMS<0)])/(nsim)
  
}
diff_total_MEMS_dat[3,6]<-mean(diff_total_MEMS)/change_MEMS_out_group$p_MEMS_results$summary_dat[1,1]

total_MEMS_out <- diff_total_MEMS_dat


## DIRECT EFFECT

diff_direct_MEMS_dat<-matrix(NA,nrow=3,ncol=6)
colnames(diff_direct_MEMS_dat)<-c(colnames(change_MEMS_out_group$f_MEMS_results$summary_dat),"Prop. of total MEMS explained")

r_names<-c("Direct MEMS - SAOM","Direct MEMS - TERGM","Diff. in Direct MEMS")
rownames(diff_direct_MEMS_dat)<-r_names

diff_direct_MEMS_dat[1,]<-c(change_MEMS_out_group$f_MEMS_results$summary_dat[1,],NA)
diff_direct_MEMS_dat[2,]<-c(change_MEMS_out_group_TERGM$f_MEMS_results$summary_dat[1,],NA)

diff_direct_MEMS <- change_MEMS_out_group$f_MEMS_results$mems_samples - change_MEMS_out_group_TERGM$f_MEMS_results$mems_samples
diff_direct_MEMS_dat[3,1]<-mean(diff_direct_MEMS)
diff_direct_MEMS_dat[3,2]<-sd(diff_direct_MEMS)
diff_direct_MEMS_dat[3,3]<-quantile(diff_direct_MEMS,.025,na.rm=TRUE)
diff_direct_MEMS_dat[3,4]<-quantile(diff_direct_MEMS,.975,na.rm=TRUE)
if(diff_direct_MEMS_dat[3,1]<0){
  diff_direct_MEMS_dat[3,5]<-length(diff_direct_MEMS[which(diff_direct_MEMS>0)])/(nsim)
}else{
  diff_direct_MEMS_dat[3,5]<-length(diff_direct_MEMS[which(diff_direct_MEMS<0)])/(nsim)
  
}
diff_direct_MEMS_dat[3,6]<-mean(diff_direct_MEMS)/change_MEMS_out_group$f_MEMS_results$summary_dat[1,1]

direct_MEMS_out <- diff_direct_MEMS_dat